TH1: Geospatial Analytics for Social Good

Take Home Exercise 1: Application of Spatial and Spatio-temporal Point Patterns Analysis to discover the geographical distribution of Armed Conflict in Myanmar

Author

Kendrick Teo

Published

September 7, 2024

Modified

September 18, 2024

TH1.1 Setting the Scene

Millions of people have their lives shattered by armed conflict every year. One of these is the Myanmar Civil War, a significant escalation of the long-running Myanmar Conflict in response to the 2021 coup d’etat.

Geospatial analtyics holds the potential to address complex problems facing society, such as this one. This study serves to discover the sptial and spatio-temporal distribution (spread) of the latest armed conflict in Myanmar by applying spatial point pattern analysis (SPPA) methods.

TH1.1.1 Loading R packages

The R packages we will use today are:

  • sf
  • tmap
  • spatstat
  • sparr
  • raster
  • maptools
pacman::p_load(sf, raster, spatstat, sparr, tmap, tidyverse)

TH1.2 The Data

Armed conflict data in Myanmar between 2010 and 2023 was downloaded from the Armed Conflict Location & Event Database (ACLED), an independent, impartial, international non-profit organization collecting data on violent conflicts and protests in all countries and territories in the world. We will be superimposing these locations with the geogrpahical boundary and subdivisions of the country, from the Myanmar Information Management Unit (MIMU).

We are interested in conflict events from 2019 onwards - starting from a year before the COVID-19 pandemic and 2 years before the 2021 coup - cumulating into the civil war of today.

TH1.2.1 Loading armed conflict data into tibble object

Here we load our aspatial armed conflict data into a tibble object.

acled_mya <- read_csv("data/2021-01-01-2024-09-16-Southeast_Asia-Myanmar.csv")
Rows: 53653 Columns: 31
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (20): event_id_cnty, event_date, disorder_type, event_type, sub_event_ty...
dbl (11): year, time_precision, inter1, inter2, interaction, iso, latitude, ...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(acled_mya)
# A tibble: 6 × 31
  event_id_cnty event_date         year time_precision disorder_type  event_type
  <chr>         <chr>             <dbl>          <dbl> <chr>          <chr>     
1 MMR66329      06 September 2024  2024              1 Political vio… Battles   
2 MMR66334      06 September 2024  2024              1 Political vio… Explosion…
3 MMR66353      06 September 2024  2024              1 Political vio… Explosion…
4 MMR66365      06 September 2024  2024              1 Political vio… Explosion…
5 MMR66368      06 September 2024  2024              1 Political vio… Battles   
6 MMR66422      06 September 2024  2024              1 Political vio… Explosion…
# ℹ 25 more variables: sub_event_type <chr>, actor1 <chr>, assoc_actor_1 <chr>,
#   inter1 <dbl>, actor2 <chr>, assoc_actor_2 <chr>, inter2 <dbl>,
#   interaction <dbl>, civilian_targeting <chr>, iso <dbl>, region <chr>,
#   country <chr>, admin1 <chr>, admin2 <chr>, admin3 <chr>, location <chr>,
#   latitude <dbl>, longitude <dbl>, geo_precision <dbl>, source <chr>,
#   source_scale <chr>, notes <chr>, fatalities <dbl>, tags <chr>,
#   timestamp <dbl>

TH1.2.2 Loading administrative boundaries into sf object

Three representations of Myanmar’s geography exist on the MIMU repository - Admin1 subdivides the country into its states and regions only, while Admin2 subdivides the country by its smaller districts, and Admin3 its townships. Further, the ACLED labels each incident with all 3 representations. For simplicity and ease of understanding, we shall use the Admin1 representation.

mmr_admin1 <- st_read(dsn = "data/mmr_polbnda_adm1_250k_mimu_1", layer = "mmr_polbnda_adm1_250k_mimu_1")
Reading layer `mmr_polbnda_adm1_250k_mimu_1' from data source 
  `/Users/kendricktty/Gits/smu_cs/is415-site/TakeHome/TakeHome1/data/mmr_polbnda_adm1_250k_mimu_1' 
  using driver `ESRI Shapefile'
Simple feature collection with 15 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
mmr_admin1
Simple feature collection with 15 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
First 10 features:
   OBJECTID          ST ST_PCODE           ST_RG          ST_MMR PCode_V
1         1  Ayeyarwady   MMR017          Region  ဧရာဝတီတိုင်းဒေသကြီး     9.4
2         2        Bago   MMR111          Region    ပဲခူးတိုင်းဒေသကြီး     9.4
3         4        Chin   MMR004           State       ချင်းပြည်နယ်     9.4
4         5      Kachin   MMR001           State       ကချင်ပြည်နယ်     9.4
5         6       Kayah   MMR002           State       ကယားပြည်နယ်     9.4
6         7       Kayin   MMR003           State        ကရင်ပြည်နယ်     9.4
7         8      Magway   MMR009          Region   မကွေးတိုင်းဒေသကြီး     9.4
8         9    Mandalay   MMR010          Region မန္တလေးတိုင်းဒေသကြီး     9.4
9        10         Mon   MMR011           State         မွန်ပြည်နယ်     9.4
10       11 Nay Pyi Taw   MMR018 Union Territory        နေပြည်တော်     9.4
                         geometry
1  MULTIPOLYGON (((95.20798 15...
2  MULTIPOLYGON (((96.17964 19...
3  MULTIPOLYGON (((93.36931 24...
4  MULTIPOLYGON (((97.59674 28...
5  MULTIPOLYGON (((97.1759 19....
6  MULTIPOLYGON (((97.81508 16...
7  MULTIPOLYGON (((94.11699 22...
8  MULTIPOLYGON (((96.14023 23...
9  MULTIPOLYGON (((97.73689 15...
10 MULTIPOLYGON (((96.32013 20...
qtm(mmr_admin1) # Plots the national boundaries of Myanmar

TH1.2.3 Creating sf data frame from aspatial data

The ACLED tibble contains coordinates, making it useful for plotting on our map as points. We can therefore use it to create an sf data frame using which we can plot our points on a map. The EPSG format of the import coordinates should be 4326, corresponding to the WGS84 Geographic Coordinate System.

acled_mya_sf <- st_as_sf(acled_mya, coords = c("longitude", "latitude"), crs = 4326)

TH1.3 A brief summary of Myanmar and the Myanmar Conflict

The second largest country in Southeast Asia at 676,578 square kilometres and an ASEAN member state, Myanmar features fertile tropical deltas in the south and the rugged Himalayan foothills to its north. It shares borders with China to the north and northeast, Laos and Thailand to the east and Southeast, and Bangladesh and India to the west and northwest. Despite having a 2800km coastline providing access to deep-sea ports and an abundance in natural resources including arable land, forests, minerals and natural gas, conflict-ridden Myanmar is considered a least developed country, with a GNP per capita of US$1144 (2011) (Myanmar Information Management Unit, n.d.) and a HDI of 0.6 (UNDP, 2024). Consequently, agriculture accounts for 36% of GDP and 60-70% of employment.

Myanmar is organised into seven states, seven regions, and one union territory - its capital, Naypyitaw, located at around the country’s geographic centre. The states - Chin, Kachin, Kayah, Kyain, Mon, Rakhine and Shan - are largely populated by the national ethnic communities, while most of the inhabitants of its regions - Ayeyarwady, Bago, Mandalay, Sagaing, Tanintharyi and largest city Yangon - are of the majority Bamar (Burmese) ethnicity. The following code chunk plots the states and regions of Myanmar into an interactive map.

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(mmr_admin1) + tm_polygons()
tmap_mode('plot')
tmap mode set to plotting

Crucially, Myanmar is one of the world’s most ethnically diverse countries, with as many as 135 different ethnic groups boasting a rich tapestry of culture and religious history (Maizland, 2022). While this may paint a rosy picture of harmony on the surface, this is far from the case in the country. While the majority Bamar enjoy a privileged position in society and have held of government and military positions, many of the country’s minority groups grapple with systemic discrimination, lack of economic opportunity and development, minimum representation in government, and abuses at the hands of the military.

The code chunk below visualises Myanmar’s ethnic breakdown in a pie chart, based on data from the US CIA World Factbook and presented by Maizland (2022).

# Plot the composition of ethnic groups in Myanmar
x <- c(2, 2, 3, 4, 5, 7, 9, 68)
y <- c("Indian", "Mon", "Chinese", "Rakhine", "Other", "Karen", "Shan", "Bamar")

y <- paste(y, x)
y <- paste(y, "%", sep="")

pie(x, labels = y, main = "Composition of Ethnic Groups in Myanmar (Maizland, 2022)", col = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7"))

# legend("topright", y, cex = 0.8, fill = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7"))

Because of this disparity and partially as a result of the divide-and-conquer strategy used by the British during colonial rule, lengthy armed conflicts have ensued since Myanmar’s independence between the ethnic groups (comprising ethnic armed organisations and smaller militias) and the Tatmadaw, the military, producing what has become known as the world’s longest continuing civil war (Rieffel, 2019), over issues from greater autonomy to control of natural resources. This did not stop even after the country underwent a democratic reform in the early 2010s - since 2017, the Tatmadaw had been mounting a brutal campaign against the ethnic Rohingya (Albert & Maizland, 2020). While the administration of then-state chancellor Aung San Suu Kyi defended their actions, restricting press freedoms in the process, rights groups and UN officials claim the Tatmadaw’s actions amounted to genocide.

TH1.4 Data wrangling

TH1.4.1 Coordinate system conversion

The EPSG area code for Myanmar is 4239, while the WGS 84-compatible code is 32647. Since we’re creating an owin object later, we need to first convert our coordinate system.

acled_mya_sf <- acled_mya_sf %>% st_transform(crs = 32647)
mmr_admin1 <- mmr_admin1 %>% st_transform(crs = 32647)
tm_shape(mmr_admin1) + tm_polygons() + tm_shape(acled_mya_sf) + tm_dots(size = 0.05)

TH1.4.2 Feature creation

Since the data is timestamped, we are able to plot and compare the frequency of conflict events over time. Before we do, though, we need to standardise the date stamps, then compartmentalise them into their year, month and day components, so that we can perform spatial-temporal point patterns analysis (STPPA) later. We will also aim to segregate events by whether or not they occurred before the coup on February 1, 2021.

A visualisation of the spatial-temporial distribution of points can be found in Annex A. Notice the consistent appearance throughout the entire study period of a large cluster of points towards the west of the country - roughly the location of Rahkine state on the border with Bangladesh.

Additionally, the ACLED dataset splits the conflict incidents in Shan state into North, South and East, and incidents in the Bago region into East and West. We need to combine these into

acled_mya_sf <- acled_mya_sf %>%
    mutate(event_date = dmy(event_date)) %>%
    mutate(DayOfYear = yday(event_date)) %>%
    mutate(Month_num = format(event_date, "%Y-%m")) %>%  # Year-Month format
    mutate(Month_fac = factor(format(event_date, "%B %Y"))) %>% # Full month
    mutate(Quarter_num = ((year(event_date) - 2021) * 4) + quarter(event_date)) %>% # Numbers quarters for STKDE
    mutate(Quarter = paste(year(event_date), "-", quarter(event_date), "Q"))  %>% # Quarter format
    mutate(admin1 = case_when(
      admin1 %in% c("Shan-North", "Shan-South", "Shan-East") ~ "Shan",
      admin1 %in% c("Bago-East", "Bago-West") ~ "Bago",
      TRUE ~ admin1
    ))
 [1] "Kachin"      "Shan"        "Magway"      "Sagaing"     "Chin"       
 [6] "Rakhine"     "Mandalay"    "Tanintharyi" "Kayah"       "Yangon"     
[11] "Kayin"       "Mon"         "Ayeyarwady"  "Bago"        "Nay Pyi Taw"
 [1] 15 14 13 12 11 10  9  8  7  6  5  4  3  2  1
[1] 15
 [1] "2024 - 3 Q" "2024 - 2 Q" "2024 - 1 Q" "2023 - 4 Q" "2023 - 3 Q"
 [6] "2023 - 2 Q" "2023 - 1 Q" "2022 - 4 Q" "2022 - 3 Q" "2022 - 2 Q"
[11] "2022 - 1 Q" "2021 - 4 Q" "2021 - 3 Q" "2021 - 2 Q" "2021 - 1 Q"
tm_shape(mmr_admin1) + tm_polygons() + tm_shape(acled_mya_sf) + tm_dots(size = 0.05)

tm_shape(mmr_admin1) + tm_polygons() + tm_shape(acled_mya_sf) + tm_dots(size = 0.05) + tm_facets(by="Quarter_num", free.coords=FALSE, drop.units = TRUE)

tail(acled_mya_sf)
Simple feature collection with 6 features and 34 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 297846.6 ymin: 1920825 xmax: 389155.2 ymax: 2653913
Projected CRS: WGS 84 / UTM zone 47N
# A tibble: 6 × 35
  event_id_cnty event_date  year time_precision disorder_type      event_type   
  <chr>         <date>     <dbl>          <dbl> <chr>              <chr>        
1 MMR11022      2021-01-02  2021              1 Political violence Explosions/R…
2 MMR10936      2021-01-02  2021              1 Political violence Battles      
3 MMR10933      2021-01-01  2021              1 Political violence Battles      
4 MMR10920      2021-01-01  2021              1 Political violence Violence aga…
5 MMR10906      2021-01-01  2021              1 Political violence Battles      
6 MMR10932      2021-01-01  2021              1 Political violence Battles      
# ℹ 29 more variables: sub_event_type <chr>, actor1 <chr>, assoc_actor_1 <chr>,
#   inter1 <dbl>, actor2 <chr>, assoc_actor_2 <chr>, inter2 <dbl>,
#   interaction <dbl>, civilian_targeting <chr>, iso <dbl>, region <chr>,
#   country <chr>, admin1 <chr>, admin2 <chr>, admin3 <chr>, location <chr>,
#   geo_precision <dbl>, source <chr>, source_scale <chr>, notes <chr>,
#   fatalities <dbl>, tags <chr>, timestamp <dbl>, geometry <POINT [m]>,
#   DayOfYear <dbl>, Month_num <chr>, Month_fac <fct>, Quarter_num <dbl>, …

TH1.4.3 Using owin to confine study area within Myanmar borders

Finally, we need to create an owin object that will be combined with any ppp objects we create later for SPPA.

mmr_admin1_owin <- as.owin(mmr_admin1)
plot(mmr_admin1_owin)

# summary(mmr_admin1_owin)

TH1.5 Preliminary analysis

TH1.5.1 Exploratory Data Analysis (EDA)

We can quickly gain some meaningful insights from our dataset using basic EDA techniques. For instance, with the following code chunk, we can retrieve the total number of fatalities from the armed violence during the study period on a quarterly basis, and print its corresponding summary statistics.

TH1.5.1.1 Getting quarterly distribution of fatalities and summary statistics

acled_mya_fatalities_by_quarter <- acled_mya_sf %>%
  group_by(Quarter_num) %>%
  summarise(`All fatalities` = sum(fatalities)) %>%
  ungroup()
acled_mya_fatalities_by_quarter
Simple feature collection with 15 features and 2 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: -208804.4 ymin: 1103500 xmax: 640934.5 ymax: 3042960
Projected CRS: WGS 84 / UTM zone 47N
# A tibble: 15 × 3
   Quarter_num `All fatalities`                                         geometry
         <dbl>            <dbl>                                 <MULTIPOINT [m]>
 1           1              803 ((-203795.3 2372607), (-201788.2 2366580), (-19…
 2           2             1939 ((-191409.1 2317222), (-174344.6 2321971), (-17…
 3           3             2698 ((-193181.1 2333786), (-191409.1 2317222), (-17…
 4           4             5542 ((-200024.3 2345207), (-199243.8 2366471), (-19…
 5           5             5301 ((-204784 2358873), (-195914.6 2357792), (-1935…
 6           6             5862 ((-206931.7 2359490), (-191409.1 2317222), (-17…
 7           7             4803 ((-206196.6 2354710), (-205672.3 2353283), (-20…
 8           8             3795 ((-206531.5 2357132), (-202593.1 2350220), (-20…
 9           9             3951 ((-199243.8 2366471), (-191409.1 2317222), (-17…
10          10             3822 ((-191261.5 2351144), (-187696.5 2324700), (-17…
11          11             3485 ((-197883.4 2339678), (-196507.9 2341782), (-19…
12          12             4792 ((-206931.7 2359490), (-204784 2358873), (-2029…
13          13             3955 ((-207135 2358896), (-206931.7 2359490), (-1978…
14          14             3982 ((-208804.4 2357274), (-201925.4 2354103), (-19…
15          15             4299 ((-197883.4 2339678), (-193173 2319902), (-1929…
summary(acled_mya_fatalities_by_quarter)
  Quarter_num   All fatalities          geometry 
 Min.   : 1.0   Min.   : 803   MULTIPOINT   :15  
 1st Qu.: 4.5   1st Qu.:3640   epsg:32647   : 0  
 Median : 8.0   Median :3955   +proj=utm ...: 0  
 Mean   : 8.0   Mean   :3935                     
 3rd Qu.:11.5   3rd Qu.:4798                     
 Max.   :15.0   Max.   :5862                     

We can establish that 3935 people die from the armed violence each quarter.

TH1.5.1.2 Plotting quarterly distribution of fatalities as a line chart

ggplot(acled_mya_fatalities_by_quarter, aes(x=Quarter_num, y=`All fatalities`)) +     
  # geom_bar(stat = "identity", fill = "red4") +
  geom_line(color = "red4") + geom_point() +
    labs(title = "Total Fatalities from the Armed Conflict in Myanmar by Quarter", subtitle = "1 = Q1 2021, 15 = Q3 2024",
         x = "Quarter",
         y = "All Fatalities") +
    theme_minimal()

There were fewer deaths from armed conflict in the first three quarters of 2021 - at the very start of the conflict - than Q4 2021 onwards.

TH1.5.1.3 Plotting proportion of fatalities by type of conflict

We can also plot the proportion of fatalities by the type of conflict. The ACLED classifies armed conflict events by mode and intensity of attack and by whether civilians are targeted. From the graph and pie chart below, the largest proportion of fatalities comes from battles, closely followed by explosions and remote violence. However, violence against civilians makes up a sizeable 11.4% of fatalities over the study period.

unique(acled_mya_sf$event_type)
[1] "Battles"                    "Explosions/Remote violence"
[3] "Strategic developments"     "Protests"                  
[5] "Violence against civilians" "Riots"                     
acled_mya_fatalities_by_event_type <- acled_mya_sf %>%
  group_by(event_type) %>%
  summarise(`All fatalities` = sum(fatalities)) %>%
  ungroup()
acled_mya_fatalities_by_event_type
Simple feature collection with 6 features and 2 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: -208804.4 ymin: 1103500 xmax: 640934.5 ymax: 3042960
Projected CRS: WGS 84 / UTM zone 47N
# A tibble: 6 × 3
  event_type                 `All fatalities`                           geometry
  <chr>                                 <dbl>                   <MULTIPOINT [m]>
1 Battles                               39059 ((-207135 2358896), (-206931.7 23…
2 Explosions/Remote violence            12473 ((-208804.4 2357274), (-206931.7 …
3 Protests                                381 ((-191409.1 2317222), (-174344.6 …
4 Riots                                   258 ((-21578.24 2371638), (275.0944 2…
5 Strategic developments                  127 ((-206931.7 2359490), (-201788.2 …
6 Violence against civilians             6731 ((-206931.7 2359490), (-204784 23…
ggplot(acled_mya_fatalities_by_event_type, aes(x=event_type, y=`All fatalities`)) +
  geom_bar(stat = "identity", fill = "red4") +
      labs(title = "Fatalities by event type in Myanmar",
         x = "Quarter",
         y = "All Fatalities") +
    theme_minimal()

total_fatalities <- sum(acled_mya_fatalities_by_event_type$`All fatalities`)
percentages <- 100 * acled_mya_fatalities_by_event_type$`All fatalities` / total_fatalities

# Set plotting area size
par(mar = c(1, 1, 2, 1))  # Adjust margins if needed

pie(acled_mya_fatalities_by_event_type$`All fatalities`, labels = paste0(round(percentages, 1), "%"), 
    col = rainbow(length(acled_mya_fatalities_by_event_type$event_type)), 
    main = "Fatalities by event type in Myanmar",
    cex = 1  
)
    
legend("topright", paste(acled_mya_fatalities_by_event_type$event_type, ":", round(percentages, 1), "%"), fill = rainbow(length(acled_mya_fatalities_by_event_type$event_type)), cex = 0.7)

We can further narrow down the scope of this analysis by limiting our data to a certain year, and to events where the civilian_targeting column is filled. The year 2022 saw the highest fatality numbers across the entire study period, so we will take a closer look.

acled_mya_fatalities_by_event_type_2022 <- acled_mya_sf %>%
  filter(year == "2022") %>%
  filter(civilian_targeting == "Civilian targeting") %>%
  group_by(event_type) %>%
  summarise(`All fatalities` = sum(fatalities)) %>%
  ungroup()

ggplot(acled_mya_fatalities_by_event_type_2022, aes(x=event_type, y=`All fatalities`)) +
  geom_bar(stat = "identity", fill = "red4") +
      labs(title = "Fatalities by event type in 2022 with civilian_targeting filled",
         x = "Quarter",
         y = "All Fatalities") +
    theme_minimal()

TH1.5.1.4 Mapping fatalities and violence against civilians

Finally, we can plot two choropleth maps displaying 2 crucial pieces of information: the locations with the greatest number of fatalities, and the locations where civilians are targeted in the violence the most in 2022. In the code chunks we will write to create the plots, we have to refer to the original tibble objects instead of the sf objects.

acled_mya_fatalities_by_admin1_2022 <- acled_mya %>%
  filter(year == 2022) %>%
  group_by(admin1) %>%
  summarise(`All_Fatalities` = sum(fatalities)) %>%
  ungroup()
combined <- left_join(mmr_admin1, acled_mya_fatalities_by_admin1_2022, by = c("ST" = "admin1"))
tm_shape(combined) + 
  tm_fill("All_Fatalities", n = 8, style = "kmeans") + 
  tm_borders(alpha = 0.5) + 
  tm_layout(title = "admin1 by number of fatalities, 2022", legend.outside = TRUE)

The above map shows Sagaing region as the location with the highest number of fatalities in 2022.

combined
Simple feature collection with 15 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -210008.6 ymin: 1072026 xmax: 724647.6 ymax: 3158467
Projected CRS: WGS 84 / UTM zone 47N
First 10 features:
   OBJECTID          ST ST_PCODE           ST_RG          ST_MMR PCode_V
1         1  Ayeyarwady   MMR017          Region  ဧရာဝတီတိုင်းဒေသကြီး     9.4
2         2        Bago   MMR111          Region    ပဲခူးတိုင်းဒေသကြီး     9.4
3         4        Chin   MMR004           State       ချင်းပြည်နယ်     9.4
4         5      Kachin   MMR001           State       ကချင်ပြည်နယ်     9.4
5         6       Kayah   MMR002           State       ကယားပြည်နယ်     9.4
6         7       Kayin   MMR003           State        ကရင်ပြည်နယ်     9.4
7         8      Magway   MMR009          Region   မကွေးတိုင်းဒေသကြီး     9.4
8         9    Mandalay   MMR010          Region မန္တလေးတိုင်းဒေသကြီး     9.4
9        10         Mon   MMR011           State         မွန်ပြည်နယ်     9.4
10       11 Nay Pyi Taw   MMR018 Union Territory        နေပြည်တော်     9.4
   All_Fatalities                       geometry
1              57 MULTIPOLYGON (((93411.72 17...
2              NA MULTIPOLYGON (((203949.9 21...
3            1158 MULTIPOLYGON (((-72918.03 2...
4             756 MULTIPOLYGON (((362696.3 31...
5            1008 MULTIPOLYGON (((309155.7 22...
6            1021 MULTIPOLYGON (((373551.3 18...
7            2488 MULTIPOLYGON (((-1717.607 2...
8            1068 MULTIPOLYGON (((208184.3 26...
9             711 MULTIPOLYGON (((364238.4 16...
10             13 MULTIPOLYGON (((220155 2248...
acled_mya_civiliansTargeted_by_admin1_2022 <- acled_mya %>%
  filter(year == 2022) %>%
  group_by(admin1) %>%
  summarise(
    count = n(),
    Civilians_targeted = sum(civilian_targeting == "Civilian targeting", na.rm = TRUE)
  )%>%
  ungroup() %>%
  mutate(`%_with_Civilian_Targeting` = (`Civilians_targeted` / `count`) * 100)
acled_mya_civiliansTargeted_by_admin1_2022
# A tibble: 18 × 4
   admin1      count Civilians_targeted `%_with_Civilian_Targeting`
   <chr>       <int>              <int>                       <dbl>
 1 Ayeyarwady    286                 44                       15.4 
 2 Bago-East     206                 47                       22.8 
 3 Bago-West     147                 31                       21.1 
 4 Chin          559                 83                       14.8 
 5 Kachin        589                126                       21.4 
 6 Kayah         542                104                       19.2 
 7 Kayin         631                114                       18.1 
 8 Magway       1620                244                       15.1 
 9 Mandalay     1315                264                       20.1 
10 Mon           523                 99                       18.9 
11 Nay Pyi Taw    68                  6                        8.82
12 Rakhine       643                165                       25.7 
13 Sagaing      5725                775                       13.5 
14 Shan-East      17                  6                       35.3 
15 Shan-North    494                130                       26.3 
16 Shan-South    439                118                       26.9 
17 Tanintharyi   958                185                       19.3 
18 Yangon       1192                256                       21.5 
combined <- left_join(mmr_admin1, acled_mya_civiliansTargeted_by_admin1_2022, by = c("ST" = "admin1"))
combined
Simple feature collection with 15 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -210008.6 ymin: 1072026 xmax: 724647.6 ymax: 3158467
Projected CRS: WGS 84 / UTM zone 47N
First 10 features:
   OBJECTID          ST ST_PCODE           ST_RG          ST_MMR PCode_V count
1         1  Ayeyarwady   MMR017          Region  ဧရာဝတီတိုင်းဒေသကြီး     9.4   286
2         2        Bago   MMR111          Region    ပဲခူးတိုင်းဒေသကြီး     9.4    NA
3         4        Chin   MMR004           State       ချင်းပြည်နယ်     9.4   559
4         5      Kachin   MMR001           State       ကချင်ပြည်နယ်     9.4   589
5         6       Kayah   MMR002           State       ကယားပြည်နယ်     9.4   542
6         7       Kayin   MMR003           State        ကရင်ပြည်နယ်     9.4   631
7         8      Magway   MMR009          Region   မကွေးတိုင်းဒေသကြီး     9.4  1620
8         9    Mandalay   MMR010          Region မန္တလေးတိုင်းဒေသကြီး     9.4  1315
9        10         Mon   MMR011           State         မွန်ပြည်နယ်     9.4   523
10       11 Nay Pyi Taw   MMR018 Union Territory        နေပြည်တော်     9.4    68
   Civilians_targeted %_with_Civilian_Targeting                       geometry
1                  44                 15.384615 MULTIPOLYGON (((93411.72 17...
2                  NA                        NA MULTIPOLYGON (((203949.9 21...
3                  83                 14.847943 MULTIPOLYGON (((-72918.03 2...
4                 126                 21.392190 MULTIPOLYGON (((362696.3 31...
5                 104                 19.188192 MULTIPOLYGON (((309155.7 22...
6                 114                 18.066561 MULTIPOLYGON (((373551.3 18...
7                 244                 15.061728 MULTIPOLYGON (((-1717.607 2...
8                 264                 20.076046 MULTIPOLYGON (((208184.3 26...
9                  99                 18.929254 MULTIPOLYGON (((364238.4 16...
10                  6                  8.823529 MULTIPOLYGON (((220155 2248...
tm_shape(combined) + 
  tm_fill("%_with_Civilian_Targeting", n = 8, style = "kmeans") + 
  tm_borders(alpha = 0.5) + 
  tm_layout(title = "admin1 by % of incidents targeting civilians, 2022", legend.outside = TRUE)

The above map shows that the highest proportion of armed conflicts in the westernmost Rakhine state involves violence against civilians.

TH1.5.2 Overarching Spatial-Temporal KDE

Next, we can try spatial-temporal kernel density estimation (STKDE) on the entire study area. As always, we start by creating ppp objects out of our sf objects, and checking for duplicates.

acled_mya_sf_quarters <- acled_mya_sf %>% select(Quarter_num)
acled_mya_quarters_ppp <- as.ppp(acled_mya_sf_quarters)
any(duplicated(acled_mya_quarters_ppp))
[1] TRUE

Since there are duplicates, we need to handle them using jittering.

acled_mya_quarters_ppp <- rjitter(acled_mya_quarters_ppp, retry=TRUE, nsim=1, drop=TRUE)
any(duplicated(acled_mya_quarters_ppp))
[1] FALSE

Finally, we join the owin object mapping the borders of Myanmar to the ppp object.

acled_mya_quarters_ppp <- acled_mya_quarters_ppp[mmr_admin1_owin]
plot(acled_mya_quarters_ppp)

We can now perform spatial-temporal KDE (STKDE).

st_kde <- spattemp.density(acled_mya_quarters_ppp)
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
summary(st_kde)
Spatiotemporal Kernel Density Estimate

Bandwidths
  h = 42816.12 (spatial)
  lambda = 0.0051 (temporal)

No. of observations
  53544 

Spatial bound
  Type: polygonal
  2D enclosure: [-210008.6, 724647.6] x [1072026, 3158467]

Temporal bound
  [1, 15]

Evaluation
  128 x 128 x 15 trivariate lattice
  Density range: [9.732081e-26, 2.087321e-10]

After performing STKDE, we are able to plot the quarterly STKDE objects between 2021 and September 2024. The full set of plots can be found in Annex A, but the plot for Q1 2021 (id 15) is of particular interest because of the larger than usual densities around the Saigaing and Chin areas, as well as near Yangon.

i = 15
plot(st_kde, quarter_numbers[i], override.par=FALSE, fix.range=TRUE, main=paste("Spatial-temporal KDE for", quarters[i]))

TH1.5.3 STKDE on Data Subset

Knowing that a significant part of the armed violence in Myanmar is due to systemic discrimination on the part of the government or military, it would make sense to repeat the process and plot STKDEs on the subset of the data concerning violence against civilians.

unique_type <- unique(acled_mya_sf$event_type)
unique_type
[1] "Battles"                    "Explosions/Remote violence"
[3] "Strategic developments"     "Protests"                  
[5] "Violence against civilians" "Riots"                     
acled_sf_againstCivilians <- acled_mya_sf %>% filter(event_type == "Violence against civilians")
acled_sf_againstCivilians
Simple feature collection with 6398 features and 34 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -206931.7 ymin: 1103500 xmax: 640934.5 ymax: 3042960
Projected CRS: WGS 84 / UTM zone 47N
# A tibble: 6,398 × 35
   event_id_cnty event_date  year time_precision disorder_type      event_type  
 * <chr>         <date>     <dbl>          <dbl> <chr>              <chr>       
 1 MMR66362      2024-09-03  2024              1 Political violence Violence ag…
 2 MMR66477      2024-09-03  2024              1 Political violence Violence ag…
 3 MMR66468      2024-09-02  2024              1 Political violence Violence ag…
 4 MMR66487      2024-09-02  2024              1 Political violence Violence ag…
 5 MMR66475      2024-09-01  2024              1 Political violence Violence ag…
 6 MMR66454      2024-09-01  2024              1 Political violence Violence ag…
 7 MMR66447      2024-08-31  2024              1 Political violence Violence ag…
 8 MMR66388      2024-08-31  2024              1 Political violence Violence ag…
 9 MMR66485      2024-08-31  2024              1 Political violence Violence ag…
10 MMR66505      2024-08-30  2024              1 Political violence Violence ag…
# ℹ 6,388 more rows
# ℹ 29 more variables: sub_event_type <chr>, actor1 <chr>, assoc_actor_1 <chr>,
#   inter1 <dbl>, actor2 <chr>, assoc_actor_2 <chr>, inter2 <dbl>,
#   interaction <dbl>, civilian_targeting <chr>, iso <dbl>, region <chr>,
#   country <chr>, admin1 <chr>, admin2 <chr>, admin3 <chr>, location <chr>,
#   geo_precision <dbl>, source <chr>, source_scale <chr>, notes <chr>,
#   fatalities <dbl>, tags <chr>, timestamp <dbl>, geometry <POINT [m]>, …
acled_sf_againstCivilians <- acled_sf_againstCivilians %>% select(Quarter_num)
acled_mya_againstCivilians_ppp <- as.ppp(acled_sf_againstCivilians)
any(duplicated(acled_mya_againstCivilians_ppp))
[1] TRUE
acled_mya_againstCivilians_ppp <- rjitter(acled_mya_againstCivilians_ppp, retry=TRUE, nsim=1, drop=TRUE)
any(duplicated(acled_mya_againstCivilians_ppp))
[1] FALSE
acled_mya_againstCivilians_ppp <- acled_mya_againstCivilians_ppp[mmr_admin1_owin]
plot(acled_mya_againstCivilians_ppp)

st_kde_againstCivilians <- spattemp.density(acled_mya_quarters_ppp)
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
summary(st_kde_againstCivilians)
Spatiotemporal Kernel Density Estimate

Bandwidths
  h = 42816.12 (spatial)
  lambda = 0.0051 (temporal)

No. of observations
  53544 

Spatial bound
  Type: polygonal
  2D enclosure: [-210008.6, 724647.6] x [1072026, 3158467]

Temporal bound
  [1, 15]

Evaluation
  128 x 128 x 15 trivariate lattice
  Density range: [9.732081e-26, 2.087321e-10]
for (i in 1: length(quarter_numbers)) {
  plot(st_kde_againstCivilians, quarter_numbers[i], override.par=FALSE, fix.range=TRUE, main=paste("Spatial-temporal KDE for violence against civilians events in quarter", quarters[i]))
}

TH1.6 Zooming in

The preliminary KDE plot identified large clusters of armed conflict in the Yangon area in the south off the Andaman sea, as well as the Saigaing and Chin areas to the central-west. These indications make good candidates for further investigation.

Additionally, we will be exploring the Rakhine state area bordering Bangladesh, and the Kra Isthmus, the narrow isthmus linking central Thailand to Peninsula Malaysia - which the overall STKDE failed to pick up. Rakhine state is the epicentre of the Rohingya crisis facing the eponymous Muslim ethnic group, which has actually been facing discrimination since the late 1970s

TH1.6.1 Fixed-bandwidth KDE

TH1.6.2 Adaptive-bandwidth KDE

TH1.6.3 Nearest neighbour analysis

TH1.6.4 Second order spatial point patterns analysis

TH1.6.4.1 F-function

TH1.6.4.2 G-function

TH1.6.4.3 K-function

TH1.6.4.4 L-function

TH1.7 Rakhine State

TH1.7.1 Fixed-bandwidth KDE

TH1.7.2 Adaptive-bandwidth KDE

TH1.7.3 Nearest neighbour analysis

TH1.7.4 Second order spatial point patterns analysis

TH1.7.4.1 F-function

TH1.7.4.2 G-function

TH1.7.4.3 K-function

TH1.7.4.4 L-function

References

  1. Albert, E., & Maizland, L. (2020, January 23). What forces are fueling Myanmar’s Rohingya crisis? Council on Foreign Relations. https://www.cfr.org/backgrounder/rohingya-crisis
  2. Maizland, L. (2022, January 31). Myanmar’s troubled history: Coups, military rule, and ethnic conflict. Council on Foreign Relations. https://www.cfr.org/backgrounder/myanmar-history-coup-military-rule-ethnic-conflict-rohingya
  3. Myanmar Information Management Unit. (n.d.). Country Overview. Retrieved 16 September 2024, from https://www.themimu.info/country-overview
  4. Raleigh, C., Kishi, R. & Linke, A. Political instability patterns are obscured by conflict dataset scope conditions, sources, and coding choices. Humanit Soc Sci Commun 10, 74 (2023). https://doi.org/10.1057/s41599-023-01559-4
  5. Rieffel, L. (2019, December 6). Peace and war in Myanmar. Brookings. https://www.brookings.edu/articles/peace-and-war-in-myanmar/
  6. Sullivan, D. P. (2021, October 21). Dire consequences: Addressing the humanitarian fallout from myanmar’s coup. Refugees International; Refugees International. https://www.refugeesinternational.org/reports-briefs/dire-consequences-addressing-the-humanitarian-fallout-from-myanmars-coup/
  7. UNDP (Ed.). (2024). Breaking the gridlock: Reimagining cooperation in a polarized world. United Nations Development Programme (UNDP).

Annex A: Full Quarterly STKDE Plots for Jan 2021- Sep 2024

for (i in 1: length(quarter_numbers)) {
  plot(st_kde, quarter_numbers[i], override.par=FALSE, fix.range=TRUE, main=paste("Spatial-temporal KDE for", quarters[i]))
}

#— OLD —

TH1.5 First Order SPPA - KDE

Kernel density estimation (KDE) serves to compute the intensity of a point distribution. It has two general steps: first to compute the point intensity, followed by spatial interpolation using a kernel function. Using KDE, we can determine the general region of Myanmar with the highest incidence of armed conflict.

TH1.5.1 Overall Fixed-Bandwidth KDE

TH1.5.2 Overall adaptive-bandwidth KDE

TH1.5.3 Quarterly fixed-bandwidth KDE

TH1.5.4 Quarterly adaptive-bandwidth KDE

TH1.6 First Order STPPA: Spatio-temporal KDE

TH1.7 First-order SPPA: Nearest Neighbour Analysis

To determine whether the plotted armed conflict events are clustered, dispersed or random, we can use the Clark-Evans test. Since tensions between Myanmar’s different ethnic groups already existed before the 2021 coup, it would be a good idea to run the test on the subsets of data both before and after the coup.

TH1.7 Second-order SPPA

Second order SPPA deals with variations in observations due to the way they interact with one another. The methods used for second order SPPA are the F-, G-, K- and L-functions. In this section, we will perform second-order SPPA on the conflict dataset, applying the F-, G-, K- and L-functions on subsets of the data before and after the 2021 coup.

TH1.7.1 F-function

TH1.7.2 G-function

TH1.7.3 K-function

TH1.7.4 L-function

Annex A: Spatial-temporal distribution of armed conflict occurrences in Myanmar between 2019 and 2023

Annex B: All quarterly fixed-bandwidth KDE plots

Annex C: All quarterly adaptive-bandwidth KDE plots